% ***********************************************************************************************
% m_commutecost.m 
% ***********************************************************************************************
clear all; 
close all;
format shortG;
% ==================== SET ESTIMATED PARAMETERS ===========================
delta=0.019; 
nu_private=0.117; nu_public=0.129;

gam_walk=0.0; 
gam_bike=0.147; gam_car=0.114; 
gam_bus=-0.307; gam_train=-0.228;

% ========= SET SHARE OF CAR OWNERSHIP 1950-75 ============================
% 1950, 1955, 1960, 1965, 1970, 1975
% Source: 
% 80 % (1987 PT) 
% 33.2% (1968 Book, 広島の都市交通の現況と将来_推計編, p191)

p_car=[0.1; 0.2; 0.3; 0.4; 0.5; 0.7]; 
Nt=length(p_car); 

% ========================= LOAD DATA =====================================
datafile = '../../data/processed data/quant/Hiroshima_DATA.csv'; 
ttfile = '../../data/processed data/quant/ODtime_allmodes_blocks1950v2_MATLAB.csv'; 
ttfile_87 = '../../data/processed data/quant/ODtime_allmodes_blocks_MATLAB.csv'; 

disp('>>> READ MAIN DATA <<<'); 
CITYDATA = readtable(datafile); 
TTDATA = readtable(ttfile); 
TTDATA_87 = readtable(ttfile_87); 

disp('>>> START SQUEEZE INFORMATION <<<')
ID = CITYDATA.geocode;   IDunique = unique(ID); 
if size(ID,1) ~= size(IDunique,1)
    disp(' >>> CHECK DUPLICATED LOCATION DATA <<<');  pause 
end 
clear IDunique 
N = size(ID,1); 

% ================== TRANSFORM TRAVEL TIME MATRIX  =======================================
ID2 = TTDATA.destination_id(1:N);
if min(ID == ID2) == 0
    disp(' >>> BE CAREFUL FOR ORDER OF LOCATIONS <<< ')
    pause
end

TT_walk=reshape(TTDATA.walktime,N,N)'; TT_bike=reshape(TTDATA.biketime,N,N)';
TT_car=reshape(TTDATA.cartime,N,N)';   TT_bus=reshape(TTDATA.bustime,N,N)';
TT_train=reshape(TTDATA.traintime,N,N)';

TT_walk_87=reshape(TTDATA_87.walktime,N,N)'; TT_bike_87=reshape(TTDATA_87.biketime,N,N)';
TT_car_87=reshape(TTDATA_87.cartime,N,N)';   TT_bus_87=reshape(TTDATA_87.bustime,N,N)';
TT_train_87=reshape(TTDATA.traintime,N,N)';


% ==================== COMPUTE AVERAGE TRAVEL COSTS =========================
disp('>>> START COMPUTING AVERAGE TRAVEL COSTS <<<')
% compute average costs in each nest (-\bar{c}^private and -\bar{c}^public) 
tcost_private=zeros(N,N);  tcost_public=zeros(N,N); 
tcost_private_87=zeros(N,N);  tcost_public_87=zeros(N,N); 
tcost_bike=zeros(N,N); tcost_bike_87=zeros(N,N); 

tcost_private = tcost_private + exp((-delta.*TT_bike(:,:)+gam_bike)./nu_private)+exp((-delta.*TT_car(:,:)+gam_car)./nu_private); 
tcost_public = tcost_public + exp((-delta.*TT_walk(:,:)+gam_walk)./nu_public)+exp((-delta.*TT_bus(:,:)+gam_bus)./nu_public)...
               +exp((-delta.*TT_train(:,:)+gam_train)./nu_public); 

tcost_private_87 = tcost_private_87 + exp((-delta.*TT_bike_87(:,:)+gam_bike)./nu_private)+exp((-delta.*TT_car_87(:,:)+gam_car)./nu_private); 
tcost_public_87 = tcost_public_87 + exp((-delta.*TT_walk_87(:,:)+gam_walk)./nu_public)+exp((-delta.*TT_bus_87(:,:)+gam_bus)./nu_public)...
               +exp((-delta.*TT_train_87(:,:)+gam_train)./nu_public); 

tcost_private = nu_private.*log(tcost_private);  tcost_public = nu_public.*log(tcost_public); 
tcost_private_87 = nu_private.*log(tcost_private_87);  tcost_public_87 = nu_public.*log(tcost_public_87); 
tcost_bike = tcost_bike + (-delta.*TT_bike(:,:)+gam_bike); 
tcost_bike_87 = tcost_bike_87 + (-delta.*TT_bike_87(:,:)+gam_bike); 


% adjustment
adj_tcost_private = tcost_private-log(0.5); adj_tcost_public = tcost_public-log(0.5); 
tcost_private = tcost_private - max(adj_tcost_private, adj_tcost_public); 
tcost_public = tcost_public - max(adj_tcost_private, adj_tcost_public); 
tcost_bike = tcost_bike - max(adj_tcost_private, adj_tcost_public); 

adj_tcost_private_87 = tcost_private_87-log(0.5); adj_tcost_public_87 = tcost_public_87-log(0.5); 
tcost_private_87 = tcost_private_87 - max(adj_tcost_private_87, adj_tcost_public_87); 
tcost_public_87 = tcost_public_87 - max(adj_tcost_private_87, adj_tcost_public_87); 
tcost_bike_87 = tcost_bike_87 - max(adj_tcost_private_87, adj_tcost_public_87);

% compute average cost if you have a car... 
tcost_car = -log(exp(tcost_private)+exp(tcost_public)); 
tcost_car_87 = -log(exp(tcost_private_87)+exp(tcost_public_87)); 

% compute average cost if you only have a bike ... 
tcost_nocar = -log(exp(tcost_bike)+exp(tcost_public)); 
tcost_nocar_87 = -log(exp(tcost_bike_87)+exp(tcost_public)); 

% ======= COMPUTE COMBINATION OF CAR AND NO-CAR ==========================
% compute average cost 
% for year 1950, 1955, 1960, 1965, 1970, 1975
tcost_mat = zeros(N,N,Nt); 
for i=1:Nt
    if i < 5
       tcost_mat(:,:,i) = p_car(i).*tcost_car + (1-p_car(i)).*tcost_nocar; 
    else
       tcost_mat(:,:,i) = p_car(i).*tcost_car_87 + (1-p_car(i)).*tcost_nocar_87; 
    end    
end

% transform average cost to commuting cost
Kappa_mat = exp(tcost_mat); 
save('commuting_cost.mat', 'Kappa_mat','TT_walk'); 
disp('>>> FINISH CONSTRUCTION KAPPA <<<')
